Josephson spectroscopy of a dilute Bose-Einstein condensate in a double-well potential 
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The dynamics of a Bose-Einstein condensate in a double- 
well potential are analysed in terms of transitions between 
energy eigenstates. By solving the time-dependent and time- 
independent Gross-Pitaevskii equation in one dimension, we 
identify tunnelling resonances associated with level crossings, 
and determine the critical velocity that characterises the res- 
onance. We test the validity of a non-linear two-state model, 
and show that for the experimentally interesting case, where 
the critical velocity is large, the influence of higher-lying states 
is important. 
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I. INTRODUCTION 

The behaviour of dilute Bose-Einstein condensates in 
multi-well potentials is remarkably rich. The tunnelling 
of a condensate through an optical lattice potential 
provides an atomic physics analogue of a Josephson junc- 
tion array. In principle, the analogue of a single Joseph- 
son junction can be realised experimentally by a conden- 
sate in a double- well potential |3|] , and this system has at- 
tracted considerable theoretical interest P~pT[. However, 
the experimental observation of the transition between 
the dc and ac Josephson effects is challenging because the 
small energy splitting associated with Josephson oscilla- 
tions means that thermal or quantum fluctuations will 
tend to destroy the effect even at the lowest achievable 
temperatures ||,0,|ll]-[l3| . While the energy splitting can 
be increased by lowering the barrier height, it then be- 
comes comparable with that of higher-lying states, and 
the accuracy of a two-mode approximation to describe 
the system PflPIJlOll becomes questionable. 

In this paper, we investigate the Josephson dynamics 
of a dilute Bose-Einstein condensate in a double- well po- 
tential, by considering transitions between energy eigen- 
states. We show that in the regime of most interest for 
experiments where the energy splittings are large, the 
influence of higher-lying states cannot be ignored. The 
paper is organised as follows. First, we calculate the 
eigenenergies of the double-well potential as a function 
of the barrier position. We show that in the vicinity 
of a level crossing the non-linearity leads to triangular 
structures in the eigenenergy curves. Similar structures 
have been shown to occur in a non-linear Landau- Zcncr 



model 1 14 1, and near the zone boundary in optical lattices 
[ p~5| , p~6[ | . Next we review the two-state model j|||[l4| , and 
show that as long as the influence of higher-lying states is 



small, the level structure can be reproduced using three 
parameters determined by numerical solution of the time- 
independent Gross-Pitaevskii equation. Finally in Sec- 
tion |V| we consider the time-dependent evolution as the 
barrier is moved, and analyse the transition between dc 
and ac tunnelling in terms of transitions between eigen- 
states. 



II. EIGENENERGIES OF THE DOUBLE- WELL 
POTENTIAL 

The eigenenergies of the double- well potential, V Xo , are 
found by numerical solution of the one-dimensional (ID) 
Gross-Pitaevskii (GP) equation 

id t m = -\v 2 m + v xo m + g\m\ 2 m • w 

Here g = Airna is the self interaction parameter, a is the 
s-wave scattering length, and n is the number of atoms 
per unit area in the y—z plane. Distance, time and energy 
are measured in terms of the standard harmonic oscillator 
units (h.o.u.), (h/muj) 1 ^ 2 , a;" 1 and Huj, respectively. In 
particular, we have concentrated on a system confined by 
the potential 
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describing a harmonic trap with a Gaussian potential 
barrier with height h and unit width. We have also car- 
ried out three-dimensional simulations, and verified that 
the essential physics is described by a ID model. This 
model facilitates an exploration of a wide range of values 
of the self-interaction, g, barrier height, h, and barrier 
velocity. 

Time-independent states of the form, %p(x,i) = 
e~ tfJ,t (p(x) , where fj, is the chemical potential, are found 
by numerical solution of the discretized system of equa- 
tions, 



fi 



(<p i+1 + <fi-i - 2(fi) /2A 2 + Vm + gtp\ 
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where ipi = ip(xi) is a real function, and A = 1/75 is the 
grid spacing. Note that in general, one must solve for 
both the real and imaginary parts of tp(xi) []l7|| . Since 
fi is also a variable we have n equations and n + 1 un- 
knowns. An additional equation is provided by fixing the 
normalisation of the wave function, 



fn 
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A solution is found iteratively by linearising equations 
(|) and (|) (with fi = (p n+1 ) 



n+l 



djj_ 

dipt 



(p) 



(5) 



where <p( p+1 ) is determined from the approximation t^ p ) 
by solving (^) using the bi-conjugate gradient method 
p8| . The iterative solution depends on the symmetry of 
the initial trial wave function Lp^\ For example, to find 
the lowest energy excited state we begin with a ground 
state solution at xq = and change the parity. Higher 
lying excited states can be found by beginning with an 
excited state solution of the harmonic trap and moving 
the barrier in from the side. 

In Fig. |l| we show the eigenenergies of the ground level 
g and the first and second excited levels, e\ and as 
a function of the barrier position xq with g = 0.5 and 
h = 12. A triangular structure appears in the upper level 
at each level crossing. For states g and e\ this structure 
is essentially the same as that discussed by Wu and Niu 
in the context of a non-linear two-state model [Q . The 
three states in the triangular region correspond to an 
anti-symmetric state with the population balanced be- 
tween the two wells, and two higher energy states with 
most of the population in either the lower or the upper 
well. The higher energy states are referred to as self- 
trapping states Any dissipative mechanism that 
damps the system towards the ground state will tend to 
destroy self-trapping, as pointed out in Refs. P,|7|Jl] 
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FIG. 1. The energy of states ground states g, and the first 
two excited states, e\ and ei as a function of the barrier po- 
sition with h — 12 and g = 0.5. In the vicinity of the level 
crossing the upper level is split into three states that form 
a triangular structure. The ground-excited state splitting is 
shown inset. 

For small xq the energy of the ground state is almost 
independent of the barrier position. This is because the 
transfer of atoms through the barrier exactly compen- 
sates for the change in potential energy. However, this 



cannot continue when all the atoms reach one side and 
the energy becomes strongly position dependent again at 
the critical displacement, x c (x c — 0.07 for the parame- 
ters in Fig. 1). 

At xq ~ 0.5 there is a second level crossing where it be- 
comes energetically favourable for the first excited state 
ei to move from the upper to the lower well. The energy 
level structure is similar to that at xq = 0, where it is 
the ground state population that moves from the upper 
to the lower well. In both cases, following the eigenen- 
ergy curve produces a large change or 'step' in the well 
populations. When considering a moving barrier these 
steps give rise to a resonance in the tunnelling current. 
In Section [y], we suggest the optimal scheme for observ- 
ing such tunnelling resonances. 

The energy splitting between the ground and lowest ex- 
cited state is extremely small for h — 12, see Fig. |l|(inset). 
From an experimental viewpoint, lower barriers resulting 
in larger energy splittings are more interesting. In Fig. |2j 
we illustrate how the triangular structures evolve as a 
function of barrier height with g — 0.5. The appearance 
of the loop structure, which coincides with the threshold 
for self-trapping, occurs at critical height h = 2.916 (see 
Section |l^) . The structure becomes more triangular as 
h increases. 
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FIG. 2. The energy of state ei as a function of the barrier 
position for g = 0.5 and barrier heights h = 2, 5, and 12. The 
energy curves for h = 2.91 and h — 2.92 are shown inset. The 
loop structure appears at h = 2.916. 

In Fig. H we show the eigenenergies for g = 0.5 and h — 
4. In this case, the loop structure appears for the first but 
not the second excited state. The appearance of a loop 
structure results in the breakdown of adiabatic following 
of an cigenenergy curve |L4|]. In Section we will see 
that this breakdown is associated with a discontinuity in 
the population difference between the two wells. 
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FIG. 3. The energy of states g, ei and ei as a function of 
the barrier position for h — 4 and g = 0.5. 

Finally in Fig. |], we illustrate the effect of increasing 
the self-interaction parameter to g = 5. In this case the 
energy splittings are a significant fraction of the harmonic 
oscillator energy level spacing and the critical displace- 
ment coincides with the position of the crossings between 
states ei and e 2 - In this regime, the influence of state e 2 
cannot be neglected, and a two-state approximation can- 
not be used to describe the system. 
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FIG. 4. The energy of states g, ei and 62 as a function of 
the barrier position for a barrier height h = 4 and g — 5. 

We stress that by solving the Gross-Pitaevskii equation 
in three dimensions we have obtained eigenenergy curves 
qualitatively similar to those shown in Figs. [I] to ^. 



III. NON-LINEAR TWO-STATE MODEL 

In this section we review the two-state model , and 
apply it to the specific case of a Bose-Einstein condensate 
in an asymmetric double- well. The wave function, ip(t), 
is written as 



xP(t) = 0i +<fo{t)ifa 



(6) 



where the zeroth-order (real) mode functions ipi, i = 1,2, 
are eigenfunctions of the time-independent GP equations 



(J) 



The mode potentials Vi and Vi are single-well potentials 
displaced to the left and the right of x = 0, respectively, 
and satisfy V\{x) = Va(— x) so that ipi(x) = ■ip 2 (—x). 
The time-dependent coefficients 4>i{t) can, in turn, be 
expressed as 



(8) 



where Ni(t) and 9i(t) are the population and phase of 
state i, with Nx(t) + N 2 (t) = 1. 

Recalling that the two mode functions satisfy 
( ipi \V Xo =o\ tpi ) = ( 1P2 \V Xo=0 \ip2 ) and redefining the zero 
of energy, we obtain the following non- linear equations of 
motion for cj>(t) = (0i(i), M t )) T '- 

i4>(t) = Hm))<t>{t) . (9) 

The non-linear Hamiltonian matrix, H(<j)(t)), is given by 



mm = \ 



-A + NE C -Ej 
-E 3 A - NE C 



(10) 



where N = N\ — N 2 is the population difference between 
the left and right sides of the well, and Ej and Eq are 
the coupling or Josephson energy and the self-interaction 
energy, respectively. The self-interaction energy of each 
mode is 



£c=.9(^l||^l| 2 |<M=S(^2||^2| 2 |^ 2 ; 



(11) 



and we define the shift in energy of the mode states due 
to the displacement of the potential barrier to be 



|A=<Vi|(V So =o-V ao )lV'i> 



(12) 



If states ipi and ip2 correspond to the population on the 
left and right of the barrier, respectively, then linearising 
yields 



x (i>i 



( dV X0 
V dx 



-axo 



(13) 



x =0 



where a < is the rate of change of the energy with 
displacement. As will be discussed below, this linear ap- 
proximation works well when the influence of higher-lying 
states is negligible. Finally, the coupling energy of the 
two- modes is 



i£j~-(v>iiHv 2 + ^ 0=0 )|</>2) 

where use has been made of the fact that 



(14) 
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dV XQ 
dxo 



life) =o 



(15) 



x =0 



The Hamiltonian, equation ([To|), is similar to that consid- 
ered by Wu and Niu in their study of non-linear Landau- 
Zener tunnelling [jlJJ . In the next section we discuss the 
eigenstates of H in relation to the numerical solutions 
discussed in Section y. 



IV. EIGENENERGIES OF THE TWO-STATE 
MODEL 

Within the two-state model, one may find the eigenen- 
ergies of the system by substituting 6i(t) — 2 (t) = et. 
The eigenenergies are given by the roots of the quartic 
equation [14| 



£ c e 3 + i(^-^j 2 -A 2 )e 2 



(16) 



The population difference between the left and right wells 
is obtained from the eigenenergies via 



N = 



A 



(E c - 2e) 



(17) 



An important feature of equation (|T^) is that there are 
four real roots when the coefficient of the quadratic term 
E c — E 2 — A 2 is positive and only two when it is neg- 
ative. For a symmetric double- well (A = 0), the addi- 
tional roots appear at the critical point where the self- 
interaction energy satisfies E c — Ej, For E c > Ej, 
the eigenenergies of the stationary states are ±Ej/2 and 
Ec/2, with the latter being doubly degenerate. The ad- 
ditional roots disappear if the displacem ent of the barrier 
is such that |A| = \2ax \ > \JE% — Ej. 

To test the applicability of the two-state model, we 
have compared the eigenenergies with those determined 
from the numerical solution of the GP equation. The 
energy curves are parametrised by three numbers: the 
splitting between the two lower levels, which is equal to 
Ey, the energy of the self-trapping states, Eq/2] and an 
energy gradient, a. In a two-state model, a is both the 
energy gradient of the self-trapping state, and that of the 
ground state for x > x c . For small g and high barriers, 
such that E c > Ej, a ~ —E c /2x c . 

Fig. |^ shows a comparison between the exact eigenen- 
ergies and the model curves for h — 4, g — 0.5. The 
values of Ec, Ej, and a are taken from the exact solu- 
tions. The energy gradient is matched to that of the self- 
trapping states at xo = (a — 3.187). This value gives 
the best agreement when comparing population dynam- 
ics (see Section |y]). However, due to the slight curvature 
of the eigenenergy curves a smaller value (a = 2.663) 
gives a better fit to the triangular structure. The agree- 
ment between the two-state model and the exact eigenen- 
ergies becomes less good as the influence of higher-lying 



states increases. For example, the model does not predict 
the almost flat position dependence of the upper levels 
in ei for h = 4, g = 5, see Fig. ||. 
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FIG. 5. Comparison between the exact (solid) and the 
two-state model (dashed) eigenenergies as a function of the 
barrier position for h — 4 and g = 0.5. 

To apply the two-state model one needs to known the 
value of Ej, Ec and a for any particular barrier height 
or non-linearity. In Fig. ^ we show how the energy split- 
tings vary with the non-linearity and the barrier height. 
The critical point, Ec = Ej, appears as either a critical 
non-linearity or a critical height depending on which pa- 
rameter is varied. Note that for large non-linearity, the 
existence of the second excited state, effectively puts 
an upper limit on the value of Ec- 
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FIG. 6. The dependence of the self-interaction energy, Ec, 
and the Josephson energy, Ej, on the self- interaction param- 
eter, g (with h = 4), and the barrier height, h (with g = 0.5). 

In Fig. 0, we show the left-right population difference 
of the self-trapping state. A population asymmetry ap- 
pears at the critical point, Ec — Ej, corresponding to 
the critical barrier height, h — 2.916 for g = 0.5. The 
asymmetry increases with increasing barrier height. As 
a function of the non-linearity, the self-trapping popula- 
tion first increases, then saturates, and finally decreases 
at large g due to the influence of e 2 . 
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FIG. 7. The population asymmetry of the self-trapping 
state as a function of the self-interaction parameter, g (with 
h = 4), and the barrier height, h (with g — 0.5). 



V. JOSEPHSON DYNAMICS 

We now apply the eigenstate picture to analyse the 
population dynamics when the barrier is moved uni- 
formly through the condensate at velocity v. The time- 
dependent GP equation is integrated using a Crank- 
Nicholson algorithm and the time evolution analysed in 
terms of transitions between eigenstates. A similar ap- 
proach was used previously to study vortex nucleation 
& 

In Fig. |8J we show both the eigenstate population dif- 
ferences and the calculated evolution of the population 
difference as the barrier is moved from the canter towards 
the right at a speed v — 4x 10~ 5 (bold line). This exam- 
ple is similar to the case studied in The eigenstate 
population differences appear as straight diagonal lines. 
For low velocities, the evolution is adiabatic, and the 
time-dependent solution follows the ground state popu- 
lation difference curve resulting in a dc tunnelling cur- 
rent. Above a critical velocity, there is a transition to a 
superposition of ground and excited states such that the 
population difference remains constant, apart from a de- 
caying oscillatory component (see inset). Subsequently, 
the excited state component encounters a level crossing 
with higher-lying states which leads to further 'steps' in 
the population difference or resonances in the tunnelling 
current. For the example of v = 4 x 10~ 5 shown in Fig. 
the speed is just below the critical velocity for the tran- 
sition from e2 to leading to a large tunnelling current 
as the object moves past xq ~ 0.9. 



A. Critical velocity 

By defining 9(t) — 9±(t) — (^(i), and rearranging equa- 
tions d^lfj|), one obtains the coupled equations 



N = Ejy/\- iV 2 sin( 
EjN 



A- NE, 



c 



cos 9 



(18) 
(19) 



y/l - N 2 

Differentiating again and substituting for N and 9 gives 



Ej{Ec^NA) 
y/l-N 2 



■ sm( 



A 2 ) 



2(1- A 2 ) 



sin 28 . 



(20) 



If N is roughly constant, this equation may be re- written 
as 9 = —deU, where U(9) describes a 'tilted-washboard' 
potential. The critical velocity is determined by setting 
the minimum gradient of U(9) equal to zero. For N w 
at t = and using A = — 2ax a = —2avt, one finds 



EjE, 



J-C'C 



2a 



3. 

4a 



(21) 



For Eq 3> Ej, one can use the linear approximation, 
a = —Eq/2x c , where x c is the critical displacement, 
which yields v c ~ Ejx c . Taking the values from Fig. [l], 
Ej = 7xl0~ 4 andx c = O.OT.gives v c - 5xl0~ 5 . Accord- 
ing to equations (|l^) and ( fl9|) the maximum population 
difference occurs at approximately v c j\[2 = 3.4 x 10~ 5 , 
which is in excellent agreement with the value determined 
by numerical integration of the GP equation. 
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FIG. 8. The population asymmetries for eigenstates g , ei 
and e-2. as a function of the barrier position with h = 12 and 
g — 0.5. The population difference as the barrier is moved 
from the centre towards the right at speeds of 4xl0~ 5 is plot- 
ted as a thick black line. An expanded view of the evolution 
for short times is shown inset. 



B. Asymmetric initial condition 

As stated above, the small critical velocity for high 
barriers (typically less than a micron per second) makes 
experimental verification challenging. Although higher 
critical velocities are obtained for lower barriers, the tran- 
sition between the dc and ac regime becomes less sharp. 
This problem can be partially circumvented by noting 
that one can induce larger population changes by start- 
ing with an asymmetric well and moving the barrier back 
through the origin. This is illustrated in Fig. || with 
a lower barrier height, h = 4. Again for low barrier 



5 



speeds the population follows the ground state distribu- 
tion, whereas for faster speeds a transition to an oscilla- 
tory current is observed. The critical velocity is a factor 
of 200 times larger than for h = 12. 
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FIG. 9. The population difference as the barrier is moved 
from Xo = 1 towards the left at a constant speed of 0.007 
(black) and 0.015 (grey) for h = 4 and g = 0.5. The corre- 
sponding population differences for eigenstates g, ei, and ei 
are indicated by the solid, dashed and dotted lines. For e2, 
the gradient of N against xo is always positive, whereas for 
ei the gradient becomes negative in the vicinity of the level 
crossing. Recalling that when the barrier is moving the tun- 
nelling current is proportional to this gradient, it follows that 
the appearance of the triangular structure in e\ is associated 
with an infinite tunnelling current and a complete breakdown 
of adiabatic evolution. 

For x < —0.6 the population difference follows that 
of the excited state, e\, therefore to observe a large pop- 
ulation difference, it is important to stop the motion be- 
fore the barrier reaches xo ~ —0.4. This restricts the 
evolution to short times or low velocities. In Fig. [l(] we 
show the population difference as a function of the barrier 
speed for both a symmetric and asymmetric starting con- 
dition. Fig. [L0](a) shows that the two-state model can no 
longer predict the correct result if the second level cross- 
ing is reached, which for t = 30 occurs when v > 0.015. 

Comparison between Fig. |l^(a) and (b) illustrates that 
a larger population difference and a better demarcation 
of the critical velocity is obtained by moving the barrier 
from the edge of the condensate inwards (asymmetric ini- 
tial condition). The numerical results indicate that the 
critical velocity for this case is similar to that for the sym- 
metric initial state. However, it is not straightforward 
to predict the critical velocity for the asymmetric initial 
condition because the relative importance of the coef- 
ficients in equation (^0|) depends on the instantaneous 
value of N(t). For h = 12, g = 0.5, the critical velocity 
is a factor of three smaller than for the symmetric initial 
state. 
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FIG. 10. Comparison of the Josephson tunnelling with (a) 
a symmetric and (b) an asymmetric initial condition. The 
population difference as a function of the barrier speed v ob- 
tained from the GP equation (solid line) and the two-state 
model (dashed line) for h = 4 and g — 0.5. In (a) the barrier 
is moved from the centre outwards for a time t — 30. For 
v > 0.015 the second level crossing is reached leading to a de- 
parture between the GP simulation and the prediction of the 
two-state model. In (b) the barrier is moved from xo = 0.95 
and stopped at xo = —0.3. The parameters Ec, Ej, and a 
for the model are taken from the GP eigenenergies. 

Finally, we consider the interesting case where the in- 
fluence of higher-lying states severely limits the appli- 
cability of the two-state model. In Fig. [ll] we show the 
population change as a barrier with height ft, = 4 is moved 
through the centre of a condensate with g = 5. For these 
parameters the critical velocity is sufficiently high that 
the level crossing to a higher-lying state is reached before 
the transition to the ac regime is completed. However, a 
large population difference for different barrier velocities 
can still be observed if the motion is halted when the 
barrier reaches xq ~ —0.5. 

The large critical velocity, v c ~ 0.08, offers the best 
potential for the experimental observation of Josephson- 
like tunnelling resonances. For example, taking a typical 
sodium condensate with a trap frequency of 20 Hz, the 
black and grey curves in Fig. [ll] correspond to a bar- 
rier formed by a blue detuned laser sheet with waist 3 
microns, moved by 5 microns in 1.5 s and 1 s, respec- 
tively. For these parameters, a model involving higher- 
lying excited states is essential to give a clear picture of 
the Josephson population dynamics. 

To restate the richness of the double-well system, it is 
worth noting that as one increases the non-linearity fur- 
ther, the transition from the ground to first excited state 
corresponds to the formation of a soliton. Consequently 
by varying only two parameters, the barrier height and 
the interaction strength, for example using a Feshbach 
resonance p0[ , one can explore the complete parameter 
space between Josephson tunnelling and soliton forma- 
tion. 
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FIG. 11. The population difference as the barrier is moved 
from position 0.9 towards the left at speed of 0.05 (black) and 
0.08 (grey) for h — 4 and g — 5. The corresponding popula- 
tion asymmetries for eigenstates g, ei, and ei are indicated 
by the solid, dashed and dotted lines. 



VI. CONCLUSIONS 

In summary, we have studied the eigenstates of a dilute 
Bose-Einstein condensates in an asymmetric double- well 
potential. We have shown that in the regime of interest 
for experiments, i.e., with a large non-linearity and a low 
barrier, a two-state description is insufficient to describe 
the tunnelling dynamics. By determining the influence 
of higher-lying states we have demonstrated that large 
population differences can be observed if the barrier mo- 
tion is halted before the upper tunnelling resonance is 
reached. 
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